{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "if not any(path.endswith('textbook') for path in sys.path):\n",
    "    sys.path.append(os.path.abspath('../../..'))\n",
    "from textbook_utils import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "csv_file = 'data/cleaned_purpleair_aqs/Full24hrdataset.csv'\n",
    "usecols = ['Date', 'ID', 'region', 'PM25FM', 'PM25cf1', 'TempC', 'RH', 'Dewpoint']\n",
    "full_df = (pd.read_csv(csv_file, usecols=usecols, parse_dates=['Date'])\n",
    "        .dropna())\n",
    "full_df.columns = ['date', 'id', 'region', 'pm25aqs', 'pm25pa', 'temp', 'rh', 'dew']\n",
    "full_df = full_df.loc[(full_df['pm25aqs'] < 50)]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(ch:pa_modeling)=\n",
    "# Creating a Model to Correct PurpleAir Measurements "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've explored the relationship between PM2.5 readings from AQS and PurpleAir sensors,\n",
    "we're ready for the final step of the analysis: creating a model that corrects PurpleAir measurements.\n",
    "Barkjohn's original analysis fits many models to the data in order to find the most appropriate one.\n",
    "In this section, we fit a simple linear model using the techniques from {numref}`Chapter %s <ch:modeling>`.\n",
    "We also briefly describe the final model Barkjohn chose for real-world use.\n",
    "Since these models use methods that we introduce later in the book,\n",
    "we won't explain the technical details very deeply here. Instead, we encourage you to revisit this section after reading {numref}`Chapter %s <ch:linear>`."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's go over our modeling goals.\n",
    "We want to create a model that predicts PM2.5 as accurately as possible.\n",
    "To do this, we build a model that adjusts PurpleAir measurements based on AQS measurements. \n",
    "We treat the AQS measurements as the true PM2.5 values because \n",
    "they are taken from carefully calibrated instruments and are\n",
    "actively used by the US government for decision making. So we have reason to\n",
    "trust the AQS PM2.5 values as being precise and close to the truth."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After we build the model that adjusts the PurpleAir measurements using AQS, we then flip the model around and use it to predict the true air quality in the future from PurpleAir measurements when we don't have a nearby AQS instrument. \n",
    "This is a *calibration* scenario.\n",
    "Since the AQS measurements are close to the truth, we fit the more variable PurpleAir measurements to them;\n",
    "this is the calibration procedure. \n",
    "Then we use the calibration curve to correct future PurpleAir measurements. \n",
    "This two-step process is encapsulated in the upcoming simple linear model and its flipped form."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we fit a line to predict a PA measurement from the ground truth, as recorded by an AQS instrument:\n",
    "\n",
    "$$ \\text{PA} \\approx b + m\\text{AQS} $$\n",
    "\n",
    "Next, we flip the line around to use a PA measurement to predict the air quality:\n",
    "\n",
    "$$ \\text{True Air Quality} \\approx -b/m + 1/m \\text{PA} $$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The scatter plot and histograms that we made during our exploratory data analysis suggest that\n",
    "the PurpleAir measurements are more variable, which supports the calibration approach.\n",
    "And we saw that the PurpleAir measurements are about twice as high as the AQS measurements,\n",
    "which suggests that $ m $ may be close to 2 and $1/m$ close to $ 1/2 $. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{sidebar} Why Two Steps?\n",
    "\n",
    "This calibration procedure might seem a bit roundabout.\n",
    "Why not fit a linear model that directly uses PurpleAir to predict AQS?\n",
    "That seems a lot easier, and we wouldn't need to invert anything.\n",
    "\n",
    "Since AQS measurements are \"true\" (or close to it), they have no error.\n",
    "Intuitively, a linear model works conditionally on the value of the variable on the x-axis, \n",
    "and minimizes the error in the $ y $ direction: $ y-(b + mx) $.\n",
    "So we place the accurate measurements on the x-axis to fit the line, called the *calibration curve*. \n",
    "Then, in the future, we invert the line to\n",
    "predict the truth from our instrument's measurement. That's why this process is also called *inverse regression*. \n",
    "Calibration is a reasonable thing to do only when the pairs of measurements are highly correlated.  \n",
    "\n",
    "As a simpler example, let's say we want to calibrate a scale.\n",
    "We could do this by placing known weights---say, 1 kg, 5 kg,\n",
    "and 10 kg---onto the scale and seeing what the scale reports.\n",
    "We typically repeat this a few times, each time getting slightly different measurements from the scale. \n",
    "If we discover that our scale is often 10% too high ($y=1.1x$), then when we weigh something in the future, we adjust our scale's reading downward by 90% ($1/1.1 = 0.9$). Analogously, the AQS measurements are the known quantities, and our model checks\n",
    "what the PurpleAir sensor reports.\n",
    "\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's fit the model. Following the notion from {numref}`Chapter %s <ch:modeling>`, we choose a loss function and minimize the average error. Recall that a *loss function* measures how far away our model is from the actual data. We use squared loss, which in this case is $ [PA - (b+mAQS)]^2 $. And to fit the model to our data, we minimize the average squared loss over our data:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{1}{n} \\sum_{i = 1}^{n} [PA_i - (b + mAQS_i]^2\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "We use the linear modeling functionality provided by `scikit-learn` to do this. (Again, don't worry about these details for now.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "AQS, PA = full_df[['pm25aqs']], full_df['pm25pa']\n",
    "    \n",
    "model = LinearRegression().fit(AQS, PA)\n",
    "m, b = model.coef_[0], model.intercept_"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By inverting the line, we get the estimate: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True Air Quality Estimate = 1.4 + 0.53PA\n"
     ]
    }
   ],
   "source": [
    "print(f\"True Air Quality Estimate = {-b/m:.2} + {1/m:.2}PA\") "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is close to what we expected. The adjustment to PurpleAir measurements is about $ 1/2 $."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model that Barkjohn settled on incorporated the relative humidity:"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "\\text{PA} \\approx b + m_1 \\text{AQS} + m_2 \\text{RH}\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is an example of a multivariable linear regression model---it uses\n",
    "more than one variable to make predictions. We can fit it by minimizing the average squared error over the data:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{1}{n} \\sum_{i = 1}^{n} [PA_i - (b + m_1AQS_i + m_2 RH_i]^2\n",
    "\\end{aligned}\n",
    "$$\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we invert the calibration to find the prediction model using the following equation:"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "\\text{True Air Quality} &\\approx - \\frac{b}{m_1} + \\frac{1}{m_1} \\text{PA}\n",
    "    - \\frac{m_2}{m_1} \\text{RH}\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We fit this model and check the coefficients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True Air Quality Estimate = 5.7 + 0.53PA + -0.088RH\n"
     ]
    }
   ],
   "source": [
    "AQS_RH, PA = full_df[['pm25aqs', 'rh']], full_df['pm25pa']\n",
    "model_h = LinearRegression().fit(AQS_RH, PA)\n",
    "[m1, m2], b = model_h.coef_, model_h.intercept_\n",
    "    \n",
    "print(f\"True Air Quality Estimate = {-b/m:1.2} + {1/m1:.2}PA + {-m2/m1:.2}RH\") "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In {numref}`Chapters %s <ch:linear>` and {numref}`%s <ch:risk>`, we will learn how to compare these two models by examining things like the size of and patterns in prediction error. For now, we note that the model that incorporates relative humidity performs the best."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
